/*
 * Phi.cpp
 *
 *  Created on: May 24, 2011
 *      Author: marchi
 */

#include "Phi.h"

Phi & Phi::operator()(Rho & myrho, bool eFilter){

	nx=myrho.Nx();
	ny=myrho.Ny();
	nz=myrho.Nz();
	std::cout << " Illo 3" <<  nx << " " << ny << " " << nz << " " << std::endl;

	Mt=myrho.getMt() ;
	matrix & oc=Mt.getOC();
	double Volume=Mt.getVol();
	int nfx,nfy,nfz;
	nfx=(nx % 2 == 0)? nx/2: nx/2+1;
	nfy=(ny % 2 == 0)? ny/2: ny/2+1;
	nfz=(nz % 2 == 0)? nz/2: nz/2+1;

	array3<Complex> ro(nx,ny,nz);
	for(unsigned i=0;i<nx;i++)
		for(unsigned j=0;j<ny;j++)
			for(unsigned k=0;k<nz;k++){
				ro[i][j][k].real()=myrho[i][j][k];
				ro[i][j][k].imag()=0.0;
			}

	fft3d Forward3(-1,ro);
	fft3d Backward3(1,ro);
	double mw1,mw2,mw3,mw,fact;
	int ia,ja,ka;
	int nx0=static_cast<int>(nx);
	int ny0=static_cast<int>(ny);
	int nz0=static_cast<int>(nz);

	std::cout << " Illo 3" <<  nx0 << " " << ny0 << " " << nz0 << " " << std::endl;

	for(int i=0;i<nx0;i++){
		ia=(i<nfx)?i : i-nx0;
		for(int j=0;j<ny0;j++){
			ja=(j<nfy)?j : j-ny0;
			for(int k=0;k<nz0;k++){
				ka=(k<nfz)?k : k-nz0;
				mw1=oc[XX][XX]*ia+oc[XX][YY]*ja+oc[XX][ZZ]*ka;
				mw2=oc[YY][XX]*ia+oc[YY][YY]*ja+oc[YY][ZZ]*ka;
				mw3=oc[ZZ][XX]*ia+oc[ZZ][YY]*ja+oc[ZZ][ZZ]*ka;
				mw=mw1*mw1+mw2*mw2+mw3*mw3;
				if(!j && !j && !k) fact=0.0;
				else fact=1.0/(mw*Volume*M_PI);
				ro[i][j][k]=fact*ro[i][j][k];

				if(eFilter) {
					ro[i][j][k]=ro[i][j][k]*filter[i][j][k];
				}
			}
		}
	}

	Forward3.fft(ro);

	Backward3.fftNormalized(ro);

	*this=ro;
	return *this;
}

Phi & Phi::operator=(array3<double> & a){
	try{
		if(this->Size() != a.Size()) throw "Phi and array3 arrays do not have the same size ";
	}
	catch(char * s){
		std::cout << s << std::endl;
		exit(1);
	}
	for(unsigned int i=0;i<nx;i++)
		(*this)[i]=a[i];
	return *this;
}
Phi & Phi::operator=(array3<Complex> & a){
	try{
		if(this->Size() != a.Size()) throw "Phi and array3 arrays do not have the same size ";
	}
	catch(char * s){
		std::cout << s << std::endl;
		exit(1);
	}
	for(unsigned i=0;i<nx;i++)
		for(unsigned j=0;j<ny;j++)
			for(unsigned k=0;k<nz;k++){
				(*this)[i][j][k]=a[i][j][k].real();
			}
	return *this;
}

void Phi::MakeFilter(int mx,int my, int mz){
	if(!filter.Size()) {
		filter.Allocate(nx,ny,nz);
		std::cout<< filter.Nx() << std::endl;
		exit(1);
	}
	fft3d Forward3(-1,filter);

	double fact=1.0/(double) ((2*mx+1)*(2*my+1)*(2*mz+1));

	for(int i=-mx;i<=mx;i++){
		int ia=(i>=0)?i:i+nx;
		for(int j=-mx;j<=mx;j++){
			int ja=(j>=0)?j:j+ny;
			for(int k=-mx;k<=mx;k++){
				int ka=(k>=0)?k:k+nz;
				filter(ia,ja,ka)=fact;
			}
		}
	}
	Forward3.fft(filter);
	for(unsigned i=0;i<nx;i++)
		for(unsigned j=0;j<ny;j++)
			for(unsigned k=0;k<nz;k++)
				filter[i][j][k].imag()=-filter[i][j][k].imag();
}

Phi operator-(Phi & a, Phi & b){
	try{
		if(a.Size() != b.Size()) throw "Phi arrays do not have the same size ";
	}
	catch(char * s){
		std::cout << s << std::endl;
		exit(1);
	}
	Phi temp(a.nx, a.ny, a.nz);
	temp=a;
	temp-=b;
	return temp;
}
Phi operator+(Phi & a, Phi & b){
	try{
		if(a.Size() != b.Size()) throw "Phi arrays do not have the same size ";
	}
	catch(char * s){
		std::cout << s << std::endl;
		exit(1);
	}
	Phi temp(a.nx, a.ny, a.nz);
	temp=a;
	temp+=b;
	return temp;
}
